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Abstract 

A variety of complex fluids consist in soft, round objects (foams, emulsions, assemblies of copoly- 
mer micelles or of multilamellar vesicles — also known as onions). Their dense packing induces 
a slight deviation from their prefered circular or spherical shape. As a frustrated assembly of in- 
teracting bodies, such a material evolves from one conformation to another through a succession 
of discrete, topological events driven by finite external forces. As a result, the material exhibits 
a finite yield threshold. The individual objects usually evolve spontaneously (colloidal diffusion, 
object coalescence, molecular diffusion), and the material properties under low or vanishing stress 
may alter with time, a phenomenon known as aging. We neglect such effects to address the simpler 
behaviour of (uncommon) immortal fluids: we construct a minimal, fully tensorial, rheological 
model, equivalent to the (scalar) Bingham model. Importantly, the model consistently describes 
the ability of such soft materials to deform substantially in the elastic regime (be it compressible 
or not) before they undergo (incompressible) plastic creep — or viscous flow under even higher 
stresses. 

PACS numbers: 83.10.GrConstitutive relations 83.80.IzEmulsions and foams in Rheology 83.50.AxSteady 
shear flows, viscometric flow 83.85.LqNormal stress difference measurements 
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I. INTRODUCTION: ELASTICITY AND PLASTICITY IN FOAMS AND EMUL- 
SIONS 



A. Prom crystals to foams and emulsions 

Historically, descriptions of deformations in crystalline, solid materials are based on a de- 
composition in terms of elastic and plastic components. Conceptual and technical problems 
arise in this process. On the one hand, general elastic formulations use continous deforma- 
tions with respect to a reference state associated with the ordered structure of minimal 



energy. On the other hand, p 
structure, called dislocations 



asticity is related to the existence of defects in the crystalline 
which are set into motion above some local threshold 
stress. Elementary motion steps constitute discrete relaxation events, which result in a drift 
of the reference state. 

In foams and dense emulsions, the local arrangement of elementary objects (droplets, 
foam cells) can be highly disordered. In the framework of crystals, this corresponds to the 
limit of a very high concentration of dislocations. Hence, a small increment of stress gives 
generally access to a large number of relaxation events. This limits the accessible range of 
purely elastic deformations between successive elementary relaxations. 

Like in crystals, each relaxation event is associated with a topological flip: the stucture 
of the network locally changes. In foams and emulsions, such events are known as "Tl 
processes": nearest neighbour links are exchanged between two pairs of cells. This class 
of materials thus exhibits an original interplay of elasticity (geometry and continuity), and 
plasticity (topology and discreteness). 



B. Immortal vs. aging fluids 



Foams and emulsions usually undergo spontaneous evolution such as coarsening (due to 



o spo 
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coalescence or ripening) or drainage [4, 15|, |6| . Such changes may induce a few topological 
rearrangements and cause the material rheological properties to slowly evolve (7|, l8| — a 
phenomenon known as aging. If the material is subjected to a weak external stress (far 
below the plastic threshold), such rearrangements may also statistically induce some creep 



which would otherwise not occur 



9|. 



In other materials made of soft, round objects, the relevant molecular processes are 



slow and no subtantial aging is observed. Among such materials — which may be called 



'immortal' — are copolymer micelles 
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11| and some foams and emulsions formulated in 



such a way as to make the diffusion of the dispersed phase and the rate of film ruptures 
imperceptible within the experimental time-scale. 

C. A brief history of flow localization 

The structural characteristics summarized above lead to interesting non-linear mechanical 
behaviours in which a rich interplay exists between structural and mechanical responses. 
One of the most extensively studied problems concerns flow localization, studied in various 
materials, from micellar solutions to granular flows. 

In the thoroughly studied system of surfactant solutions self-organized as giant micelles, 
the flow curve was originally observed to exhibit a plateau-like behaviour under controlled 
shear rate [l^ . 

Structural observation followed, demonstrating shear-banding. Thus, in the stress plateau 
region, two different organizations of the material coexist: an isotropic region, similar in 
structure to the original solution, and a strongly birefringent region, in which the micelles 



are aligned to a high degree with the flow direction [13|. 

This situation was initially interpretated in terms of out-of-equilibrium phase transitions 
in the material, leading to a steady-state coexistence between two structurally homogeneous 
domains jl^ . 

Theoretical descriptions attempted to capture the onset of this transition and the stability 



of the coexistence in terms of interfacial dynamics and mechanical instability 



i3- 



A more detailed investigation of the birefringent phase has more recently revealed strong 



spatial [16|] and temporal [17|] variations. That disagrees with the initial simple picture. 

The birefringent phase consists in numerous transient, narrow zones of very large shear. 
The term "fluid fracture" has been proposed to describe these individual events, which 
have been observed in different systems with similar rheological properties but different 



internal structures (connected microemulsions [19|, copolymer cubic-phase jlOl. 



The understanding of the shear-banding phenomenology has thus begun to shift from a 
phase transition picture to a fracture picture, and the interest is now focusing on the highly 
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localized events that initiate the transition, and on their relation to the structural properties 
of the material. 

Similarly, in the case of foams, a similar phenomenology of localization has been observed 



both in experiments 



20| and simulation 2l|, and interpreted j4J] in terms of the interaction 



between individual events via (mostly elastic) deformations of the surrounding material 
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2J|. The challenge now consists in understanding the self-organization of dispersed 



relaxation events into a localized fracture-like behaviour. 

Thus, the emergence of a fluid fracture from these discrete, elementary relaxation events 
appears as a well-defined problem. In that respect, two main problems remain open: (i) 
the role of structural disorder, and (ii) non-local effects between individual events, mediated 
through elastic stresses in the material Both problems relate to the unknown typical length 
scale on which the discrete system should be averaged for a descriptions in terms of a 
continuous model. 

More precisely, problem (i) addresses this length scale from the limit of smaller length 
scales where disorder is relevant. Conversely, problem (ii) addresses it from the larger 
length scale limit: a suitable constitutive equation, incroporated into the general framework 
of continuum mechanics, provides the tools for generating such non-local effects. 



D. Ingredients of our model 

In the present work, we focus on the second problem discussed above, and construct an 
example of a rheological model inspired by the behaviour of such "immortal" fluids. It is 
characterized by four main features. 

1. The flow properties are motivated and discussed in terms of microscopic considerations 
(Tl processes). 

2. In order to incorporate the non-local elastic effects mentioned above, our model is 
written in a fully tensorial form, whether in two and in three dimensions. 

3. It is a commonly observed feature that such soft materials deform substantially before 
they display plasticity. In other words, their yield stress is comparable to their elastic 
modulus (unlike for classical, hard crystals). Correspondingly, the present model im- 
plements a consistent description of the elastic properties of the material up to finite 
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deformations {i.e., beyond the usual approximations valid at small deformations). 

4. As we shall see now, an amorphous, elastic material undergoing plasticity loses the 
memory of past events. The initial reference state thus has no physical relevance. 
Correspondingly, our model is developed in the Eulerian formalism (attached to the 
current reference state). 



II. LOSS OF MEMORY AND CONSEQUENCES 





FIG. 1: When a block of foam undergoes very large deformations (symbolized by the evolution 
of the coordinate system between the three left-hand side drawings), one might think that its 
local structure correspondingly evolves towards a highly stretched configuration. In fact, it always 
remains similar to its initial state, as illustrated on the right-hand side drawing. 

In foams or emulsions and in crystals alike, large deformations of the sample do not 
imply large deformations of individual objects, since topological rearrangements relax local 
stresses. In both types of systems, two objects that are initially in contact can end up 
at a large mutual distance when many topological rearrangements have occurred. In such 
a situation, there does not exist any kind of elastic restoring force between both objects. 
Hence, the distance between them is irrelevant to the current mechanical behaviour. As a 
consequence, the deformation from the initial state, which keeps track of such large distances, 
is mechanically irrelevant, even though it has a clear experimental meaning. In other words, 
the material has lost the memory of such large deformations. 

As we shall see, this is the reason why: 

• we use Eulerian coordinates; 

• we define a quantity called "stored deformation"; 
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• we discuss the impact of local disorder. 
A. Eulerian description 




FIG. 2: Coordinate systems for a deformed foam. Let us attach a system of coordinates to the 
initial state of a foam (left-hand side drawing). After some deformation, we may describe the 
quantities attached to the foam in terms of either of two coordinate systems (right-hand side 
drawings): either the initial (now deformed) coordinate system, or a new (undistorted) coordinate 
system defined on the current state of the foam. In the case of elastic deformations, when material 
keeps trace of its initial configuration, such a choice does not have noticeable consequences. But 
in the case of deformations that imply plastic events which progressively erase the memory of the 
initial state (see Figure [1]), the choice of the initial, much deformed coordinate system would not 
be physically (or computationally!) particularly relevant. 

In order to describe the deformations of a material sample and the evolution of the phys- 
ical quantities attached to it, two types of coordinate systems are commonly used: either 
"Lagrangian" coordinate systems attached to the initial state of the sample, or "Eulerian" 
systems attached to its current configuration, which coincide to first order in a small defor- 
mation context. 

In the case of an elastic material, that keeps the memory of its initial state, the choice 
of one or the other does not have consequences other than computational (physicists most 
usually use Lagrangian coordinates jl|). 

However, in the case of a material that progressively looses the memory of its initial 
configuration, such as a foam or an emulsion via rearrangements (Tl processes described 
later in this paper), it would be physically irrelevant (and computationally tedious) to refer to 
the initial sample state. One therefore generally uses Eulerian coordinates in such situations. 
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For instance, the Navier-Stokes equations are usually expressed in a fixed coordinate system 
(Eulerian approach). 



B. From deformation to 'stored deformation' 



Experimentally, the accessible variables are (1) the deformation or deformation rate, as 
measured or imposed at the sample boundary, and (2) the stress (at least in some systems 
such as foams, where it can be extracted from the shape of the individual objects, or in 
photoelastic systems). In order to set up a spatio-temporal numerical scheme, one needs 
not only continuity and force balance equations, but also a specific evolution relation of the 
local stress in terms of the deformation rate. Our purpose in this paper is to propose a 
model example for this missing ingredient, and to explore its properties. 

The very notion of a fixed reference state, and of a global deformation from this state, 
being of no use, the stress, as an index of the local elastic deformation, represents only 
the recoverable (or 'stored') part of the deformation. It can be defined through the fol- 
lowing thought experiment, described on Fig. [31 a fragment of the material is cut, in the 
deformed state, and allowed to relax; the "stored deformation" is defined as the inverse of 
the deformation observed during this relaxation. 



In the case of foams, it was shown a few years ago 



25l ] that it is possible to construct a 



deformation tensor from the experimentally observed inter-bubble (centre-to-centre) vectors. 



which indeed faithfully represents the stress 



26|. 



C. Disorder, stored deformation, and stress 

The knowledge of the local stored deformation is exactly equivalent to the knowledge of 
the stress. They are related through a specific, material-dependent, constitutive relation, 
namely the elasticity, be it linear or not. 

Generically, disordered systems are locally frustrated and contain internal stresses even 
in the absence of external applied stress. In other words, stored deformations are nonzero. 
One might think it possible to relax stored deformations by cutting the material into pieces 
and sewing them together again, as described above. In fact, the relaxed pieces do not 
fit together nicely, even after adjusting local orientations: the field of stored deformations 
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FIG. 3: Stored deformation: operative definition through a thought experiment. Consider stretched 
piece of material Q, When its macroscopic deformation (represented by a large ellipse) is relaxed 
(center, large circle), there remain local internal stresses. Only by cutting out small pieces of 
material (right) can the corresponding stored deformations (small ellipses) be relaxed (small cir- 
cles). The resulting pieces ri* cannot fit together without restoring local stresses: relaxation is 
meaningful only on the local scale. Note that the orientation of each piece is arbitrary and can 
always be chosen in such a way that when going back to the initial, macroscopically stretched state 
(left), the corresponding transformation (F, see Eg. Ilip be purely elongational (no rotation), with 
stretching factors Aj (see Eq. [T2]) . This stored deformation is related to the local stress via the 
material elasticity (see Eq. 

cannot be reconstructed from a displacement field. If one were to sew all pieces together 
again, one would need to stretch each of them appropriately, thus reconstructing a frustrated 
stress field corresponding to a state with zero external applied stress. 

The notion of reference state is always clear in a local context. But it cannot be extended 
to any macroscopic part of a disordered material. Indeed, it would not be extensive: one 
half of a relaxed sample generally does not match the relaxed state of the same sample half. 
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D. Evolution of the stored deformation 



Let us return to the construction of the local evolution of the stress in terms of the applied 
deformation rate. 

The main point is that topological events participate in the applied deformation but relax 
part of the corresponding stored deformation. The evolution of the stored deformation must 
thus involve both an entrainment part and a plastic part. 

The entrainment part is purely kinematic, driven by the velocity gradient. The plastic 



part reflects the Tl relaxation processes triggered at large stored deformations. It always [45] 
tends to lower the stored deformation. It reflects the rate at which the material looses 
memory of the local reference state which is implicit in the stored deformation. 

III. CHOICE OF A RHEOLOGICAL MODEL 

Let us now choose a rheological model. The considerations of Section fllDI the rheological 
model of a foam must incorporate a spring (which represents elasticity) in series with a 
creeping, plastic part. In order to choose this plastic part, let us now review a few common 
rheological models. They correspond to the generic form given by Figure HI where we added 
an optional viscous element in parallel with the other two elements. Note that such a viscous 
element impacts the stress response of the material, but not the local dynamics of stored 



deformation (at least not directly 461]). 

Such models are listed on Table IIIII together with an indication of their creep and relax- 
ation properties. We now review some of these models, which have been used in the context 
of foams or similar materials. 



A. Burger model for weak applied stresses 



Is 



The rheology of dry liquid foams under weak stresses is well-described j9[ by the Burger 
model, which consists in a Maxwell group in series with a Kelvin- Voigt group (see Table UTTll . 
The elastic response (resulting from both spring elements) corresponds to the deformation 
of the disordered network of inter-bubble films and Plateau borders. The (short) transient 
(given by the Kelvin- Voigt group) corresponds to the viscous stretching of films needed 
to reach the new equilibrium film conformation. The (slow) creep (given by the Maxwell 
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FIG. 4: The considerations of Section [11 Dl lead us to the idea of a rheological model consisting in a 
spring in series with some creeping element yet to be defined. The long, soft spring reflects the fact 
that creep may trigger at rather large elastic (stored) deformations. Optionnally, an additional 
viscous element may be added in parallel: in a context of imposed deformation, it will not alter 
the dynamics of the system. 

viscous element) corresponds to the spontaneous Tl processes (which are responsible for 
aging phenomena) being slighlty biaised by the ambiant stress. 

B. Models with non-zero threshold 

Although well suited to describe the response of foams at weak stresses, the Burger model 
does not incorporate the finite plasticity threshold. Let us now review some models that 
do incorporate the threshold, even though they oversimplify the short-time response under 
weak stresses. 

More precisely, these models are designed to provide most, if not all, of the following 
features: (i) a. simple elastic response to small stresses, (ii) yielding above a stress threshold, 
and (in) a viscous response at large, constant deformation rates. 

The first family of such models includes the simple elasto-plastic model (a spring in 
series with a solid friction element), with an optional viscous element in parallel (Table UTTl) . 
The simple elasto-plastic model does not provide feature (in). The viscous element was 



incorporated by Marmottant and Graner 



271] to account for this feature, observed in foams. 



The second family, also used in the context of foams or rheologically similar materials (see 

ement coupled to the solid friction element, i.e., the Bingham 
291], and present work), also with an optional viscous element 



Table IIIip , includes a viscous e 
model (Takeshi and Sekimoto 



in parallel (Saramito 281]). 



In all cases except the pure elasto-plastic model, the viscous element provides feature (Hi). 
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Model 



Relaxation 
from fj > fj 



y 



towards a = a,. 



Creep 
threshold 



Foams 
(or similar) 



* AWV"* LC2 • linear Maxwell 


delayed 


ay = 
(viscoelastic) 




D—^"^'**'^ " ^ linear Burger 


delayed 


ay = 
(viscoelastic) 


Hohler [9] 


P linear elasto-plastic 


immediate 


ay-^G 


Marmottant-Graner [27 
(+ viscous) 




linear Bingham 



delayed 



ay-^G 



Saramito [28] (+ viscous) 
Takeshi-Sekimoto [29.1 



• — vwwvvw 



non-linear Bingham delayed 



G 



present model 



TABLE I: Some common rheological models with creep. Each model is schematically designated by 
a combination of springs, viscous and frictional elements, where viscous elements are not necessarily 
assumed to respond linearly. The ability of the material to relax from above the threshold stress 
when applied deformation is stopped is indicated. The value of the stress threshold is compared 
to that of the elastic modulus in order to estimate the deformation at the threshold. Models have 
been labeled as "linear" when the deformation at the threshold is small, whether explicitely or 
implicitely. Some references are given when such models have been used in the context of foams or 
other soft, disordered materials. In some instances (marked with label "+ viscous"), an additional 
viscous element was introduced in parallel with the other elements altogether. This element is then 
indicated with dashed lines in the corresponding diagramme in the left-hand side column. 

C. Each foam has its own rheological model 

The reason why different models have been suggested is that the rheology of a foam 
varies with several parameters, among which surface tension, bubble size, polydispersity, 
surfactant properties. Let us concentrate on the effect of volume fraction and continuous 
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phase viscosity, see Figure [51 



1. Elasticity and volume fraction 

At gas volume fractions below the close-packing threshold, the foam flows and displays 
no elasticity. 

Slightly above the threshold, elasticity is weak at small stored deformations as bubbles 
move rather freely between neighbouring bubbles, and strengthens at larger stored deforma- 
tions as they come into closer contact with neighbours. 

At large gas volume fractions, the elastic modulus is expected to be large even at small 
deformations, as the bubbles are already in close contact. 



2. Plastic threshold and volume fraction 

In order for neighbouring bubbles to undergo a plastic rearrangement (such as a Tl 
process), they need to deform more importantly when the gas volume fraction is larger. As 
a result, the plastic threshold ay is also expected to be larger. 



3. Large elastic deformations 

The plasticity threshold of foams usually corresponds to moderate deformations, that 
are beyond the small deformation regime. For instance, the threshold deformation for a 
polydisperse foam under shear in the dry limit (volume fraction approaching unity) is on 



the order of 30 to 50% 



30|. 



The ratio Cy/ G is an indication of whether elastic non-linearities appear prior to the onset 
of plasticity. It is not clear to us whether this ratio increases or decreases with gas volume 
fraction. 



4- Relaxation and volume fraction 

Under stationary conditions, elasticity is inactive and both viscous elements play similar 
roles. They can be distinguished, however, in transient responses. 
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At moderate volume fractions, we expect any relative motion of bubbles to generate 
viscous dissipation: the general viscous element should dominate over the viscous friction 
element. 

At higher gas volume fractions, bubbles interact so intimately that most viscous dissi- 
pation can be expected to arise during plastic events. In other words, the viscous friction 
element should dominate over the general viscous element. 

The relative weight of the viscous friction element and the general viscous element is 
apparent on Figure [51 

Their ratio also impacts the ability of the material to relax when the applied deformation 
rate is suddenly brought to zero. Indeed, in this respect, only the viscous friction element is 



relevant. A model without such a viscous friction element 



27| does not display relaxation in 



common situations (like oscillatory measurements) where the deformation rate may reverse. 



5. Viscous vs. plastic behaviour 

Besides, the visosity of the continuous phase impacts the relative importance of the 
viscous elements and the solid friction element. This impacts the stationary response of a 
foam, which typically changes from mainly plastic to mainly viscous as the deformation rate 
is increased. This transition is expected to occur at lower deformation rates if the continuous 
phase viscosity is increased. 



D. Choice of the Bingham model 

Let us now choose a specific model in order to develop a fully tensorial version of it. 

Except for sollicitations at very low deformation rates, the viscous elements impact the 
rheological response. Should we keep them both? 

The general viscous element provides an additional contribution to the stress. It does not 
represent any technical difficulty. Besides, it does not impact the behaviour of the system 
under imposed applied deformation conditions. 

By contrast, we believe that the viscous friction element provides essential features such 
as the ability to display relaxation. 

For these reasons, in the remaining part of this paper, we focus on the Bingham model 
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FIG. 5: Tentative variation of foam rheology as a function of its volume fraction and of the 
continuous phase viscosity. Four instances of the local aspect of the foam are drawn in the center 
of the figure, for low (light grey) or high (dark grey) continuous phase viscosity and for dispersed 
(round) or concentrated (faceted) bubbles. In each case, a tentative corresponding rheological 
model is schematically represented in terms of one spring, one solid friction element, one viscous 
friction element, and one general viscous element. The strength of each element is coded as weak 
(dashed line), medium (thin solid line) or stroqg (thick line). 



(see Fig. [6]), which is the simplest one to provide all three desired properties together with 
relaxation. We hope that it may also apply to a broad range of materials made of densely 
packed, soft, essentially round objects. 

We shall keep in mind that the parameters of the model (G, (jy and 77 on Fig. [H]) will 
depend on such physical quantities as the volume fraction and the viscosity of the continuous 
phase, as discussed in Section UlI CI and illustrated by Fig. O 

E. Behaviour of the Bingham model 




FIG. 6: With the Bingham model, from which we develop a fully tensorial model in the present 
work, the response to weak stresses is elastic (modulus G); by contrast, the response to large 
stresses presents both a yield threshold (dy) and a viscous component (r?). 

Let us now check (with its scalar version) that the Bingham model (Figure E]) behaves 
as expected for this class of materials. In this model, the evolution of stress (for positive 



values of the stress 



4 Ti p is given by: 



a = GD-^^-j^e\a-ay] (1) 

where 9{x) = 1 when x > and d{x) = otherwise, and where D is the applied deformation 
rate, i.e., the symmetric part of the velocity gradient tensor: 

O.Yl^ (2) 

Equivalently, the evolution of the spring elongation e = cr/G is given by: 

e = D-^-^^0[e-ey] (3) 
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where Sy = <7y/G. 

Let us now consider successively three simple experiments: (a) quasistatic imposed de- 
formation; (b) constant, imposed deformation rate; (c) constant, imposed stress. The corre- 
sponding evolution of the main rheological variables is schematically represented on Figure [3 

1. Quasistatic imposed deformation 

Under low applied deformation, the stress depends linearly on deformation. This low 
deformation regime is valid as long as the resulting stress is smaller than the yield value ay. 
At larger deformations, the stress remains constant and equal to ay, even under arbitrary 
large (but constant) deformations. 

The reason why foams and emulsions display such a solid friction behaviour is that the 
Tl processes are very similar to the relaxation of surface bumps involved in the friction 
between rough solids: in both cases, part of the mechanical work done by the imposed stress 
is dissipated in discrete relaxation events which enable discrete deformation steps. As a 
result of these events, the work is proportional to the total deformation (rather than to the 
velocity, as in a viscous fluid). 

Under imposed deformation, supposing that the deformation value is reached through a 
quasi- st at ionnary process, the stress will remain constant as soon as the elastic stress associ- 
ated with a deformation increment is exactly compensated by the stress relaxed through the 
plastic processes. In this situation, even if the imposed deformation can be arbitrarily large, 
the stored deformation remains equal to the value that corresponds to the yield stress: any 
extra applied deformation is relaxed through Tl processes. The noisy aspect of the stress 
plateau reflects the disorder of the material — and hence, of the distribution of available 
relaxation processes. 

2. Constant imposed deformation rate 

Under constant applied deformation rate (see Figure Wp), the stress rises linearly at 
short times, as long as it is smaller than ay. At later times, it eventually stabilizes above 
the threshold, and its final value a^o increases with the deformation rate (affinely in the 
Bingham model). 
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3. Quasistatic imposed stress 



Under constant applied stress, by contrast (see Figure [7t), the system displays two dif- 
ferent behaviours. At small stress values (cr < ay), it behaves elastically. It switches to a 
flow behaviour at larger stress values (a > Cy). 

k O-oo 




£ t t 



(a) (b) (c) 

FIG. 7: Behaviour of the Bingham model (time t, deformation e, deformation rate D, stress 
cr) in three series of experiments: (a) quasistatic imposed deformation (with a zoom on the stress 
fluctuations due to individual Tl processes and corresponding elastic loading periods); (b) constant, 
imposed deformation rate (top: value of plateau stress); (c) constant, imposed stress. 

F. Two-phase fluid 

1. Evolution modes 

Since a foam (or an emulsion) is a system with two different phases, it has more deforma- 
tion modes than a monophasic fluid. Figure M depicts three isotropic modes obtained from 
the initial configuration (0) via gas diffusion (1), via an applied isotropic (2) or deviatoric 
(3) stress, or via fluid permeation (4). 

Note that the gas diffusion mode and the liquid permeation mode are plastic, even for 
small magnitudes: they are accompanied by dissipation and they lead to situations that are 
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FIG. 8: Four changes in the local structure of a foam. (0) Initial structure. (1) Structure obtained 
through diffusion of gas from neighbouring regions into the bubble: the number of gas molecules 
in the bubble (represented by black dots) has increased, while the amount of liquid that surrounds 
the bubble has not changed. (2) Structure obtained by isotropic stretching: the quantity of gas 
and the quantity of liquid have not changed. This is not a stress-free conformation. (3) Structure 
obtained by deviatoric stretching: the volume, the quantity of gas and the quantity of liquid have 
not changed. This is not a stress-free conformation. (4) Structure obtained by permeation of 
liquid from the vicinity of the bubble towards other places in the foam: the amount of liquid that 
surrounds the bubble has decreased, while the number of gas molecules has not changed. In the 
present work, only changes where both phases are transported simultaneously, such as (0) (2) 
and (0) — > (3), are considered. 
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stable even in the absence of any extra applied stress. Indeed, in situation (1), the extra 
amount of gas occupies the extra volume. Similarly, in situation (4), the amount of gas in 
the bubble and the bubble volume have remained constant, so the force balance within the 
material has not been altered despite the loss of liquid. 

In the present work, for the sake of simplicity, we restrict to changes in the applied stress 
i.e., modes (0) (2) and (0) (3). For small magnitudes, these two modes are elastic: if 
the applied stress is removed, the system returns to state (0). For larger amplitudes mode 
(0) (3) is plastic. Indeed, under large applied deviatoric stresses, rearrangements of the 
liquid films between bubbles (like the Tl processes described in Section IVI Al below) lead 
to stress-free states that differ from the initial state, even though they may locally be very 
similar, if not identical, to state (0). 



2. Density, velocity and stress 

With this choice of modes (0) (2) and (0) (3) rather than modes (0) (1) or 
(0) (4), the weight fraction of each phase remains constant, and the material can safely 
be handled like a one-phase fluid, which substantially simplifies its description. 

The overall material density, p, can be expressed in terms of the contributions from both 
phases: 

P = V'Pliq+ (1 - ^)Pgas (4) 

1.^ + 1^ (5) 

P Pliq Pgas 

where ip (resp., ip"^) is the volume fraction (resp., weight fraction) of the liquid phase. 

Because we ignore the gas diffusion and liquid permeation modes depicted on Fig. [HI the 
velocity of the material at larger length scales suffers no ambiguity, as it is the same in the 
dispersed phase (bubbles or droplets) and in the continuous (liquid) phase: 

^ (^'dispersed phase) (^continuous phase 

As for the stress, it varies strongly at the microscopic scale within such a structured 
medium as a foam: compressive within a bubble and across the gas/liquid interfaces, tensile 
along these interfaces (surface tension), tensile within the liquid Plateau borders and vertices. 
In the present work, the stress variable a represents the sum of these contributions averaged 
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(6) 



at some larger length scale, where such structural details are smoothed out: 

^gas/liqinterf. 



(7) 



G. Compressibility 

In Section IIIID[ we motivated our choice of the Bingham model to describe the rheo- 
logical behaviour of an emulsion or a foam. The above considerations, however, were based 
essentially on scalar arguments. 

One of the main properties that reflect the actual, tensorial nature of deformations and 
stresses in the material is compressibility, which is the ability of the material to adapt its 
volume when the pressure {i.e., the isotropic part of the stress) is changed (see mode (0) — >■ 
(2) in Fig. [8]). A material is considered incompressible when the elastic modulus involved 
for changes in volume is much larger than the elastic modulus involved for deformations at 
fixed volume. In practice, the deformations of such a material therefore obey a fixed volume 
condition (mode (0) — > (3) in Fig. [8]). 

Among the materials we address in the present work, emulsions can be considered incom- 
pressible for all practical purposes since both phases are liquid. 

Foams can also be considered incompressible as long as their constitutive bubbles are 
not too small. Indeed, their compression modulus is then typically equal to the pressure 
in the dispersed, gas phase, while their shear modulus (and other moduli corresponding to 
deformations at constant volume) are on the order of the interphase surface tension divided 
by the typical bubble radius. Hence, under atmospheric pressure and with usual liquids, 
gases and surfactant molecules, foams with bubbles not smaller than 0.1 mm in size can 
safely be considered incompressible. By contrast, the compressibility of foams made of 
micron-sized (or even smaller) bubbles cannot be neglected. 

As for bubble monolayers, when regarded as two-dimensional foams, their apparent com- 
pressibility depends on the boundary conditions. When such a monolayer is squeezed be- 
tween two solid plates, it can be considered incompressible under the same conditions con- 
cerning the typical bubble size as a three-dimensional foam. In other situations, a bubble 
monolayer may have at least one free interface, for instance when it floats on a bath of the 
liquid, continuous phase, and/or when its upper surface is in contact with the atmosphere. 
In such a situation, each bubble is free to slightly deform in the vertical direction in order 
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to better accomodate in-plane stresses. As a result, the bubble monolayer may appear com- 
pressible as seen from above, even though the total bubble volume in fact remains essentially 
constant. 

In the present work, in order to be able to describe the rheology of all such systems, 
we consider a compressible material (modes (0) (2) and (0) — ^ (3) both available in 
Fig. [8]). Nevertheless, specific properties or mathematical formulations suitable for the 
incompressible case (mode (0) (3) only) are provided whenever appropriate. 



IV. GENERAL FORMULATION FOR MATERIALS CAPABLE OF CREEP 



A. Evolution of the stored deformation 

After the scalar description given in Section IIII El we now turn to a tensorial version of 
the model. In particular, instead of the scalar deformation e, we will now use a deformation 
tensor e, to be defined in detail in Section |Vl 

As mentioned at the beginning of Section llllt we are interested in materials that display 
some elasticity and are capable of creep (see Figure H]). The evolution of their deformation 
can be decomposed into an elastic part and a creeping part 
a decomposition was provided above for the (scalar) Bingham model, see Eq. ([3]). More 
generally, the evolution of the stored deformation can be written in the form: 

e = kinematics(e, Vv) — creep(e, D), (8) 

where D is the applied deformation rate as defined by Eq. ([2]). 

In the Bingham model, the creep term depends only on the stored deformation e. In other 
models 27|, it additionally depends on D. This point will be discussed later, in Section |XB[ 
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B. Stress evolution 

For a material capable of creep (see Section IIV Al) . once the evolution of the stored 
deformation is known, the stress evolves as prescribed by elasticity (see Figure [9]). 

Hence, in the elastic regime, since stored deformation evolves purely kinematically, the 
stress evolution results from the kinematics and elasticity alone. In the presence of plasticity, 
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kinematics 
e ► de / dt 
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FIG. 9: Evolution of stored deformation e and stress a in the elastic regime. The evolution de/dt 
of the stored deformation e is purely kinematic, as it directly results from the convection of the 
material the velocity gradient Vv. As for the stress a, it is related to the current stored deformation 
e through the elastic law. As a result, by composition, the evolution da/dt of the stress results 
from the expression of de/dt and from the elastic law. 

although the evolution of the stored deformation includes an additional term (see Eq. [8l), 
the stress evolution is still deduced therefrom in the same way, namely via elasticity. 



C. Consistency of some commonly used stress evolutions 



It is common habit to express the stress evolution directly in the form 

da 



dt 



(9) 



At first sight, this might seem exactly equivalent to the evolution for the stored deformation 
as given by Eq. ([8]). 

In fact, if no special care is taken, such an expression does not usually correspond to the 
composition of the stored deformation evolution with elasticity as illustrated on Figure [9l 
The validity of such an evolution equation for stress is then implicitely restricted to the 



domain of small stored deformations, where elas 

As an example, as shown in another work 
the Gordon-Schowalter derivative interpolation 



;icity is linear. 



32 1 , generalized Maxwell models that involve 
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33| suffers such restrictions (except for 



the special cases of upper and lower-convected derivatives). 
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V. RELAXED, REFERENCE STATE, AND ELASTICITY 




FIG. 10: Variation of the local material stored deformation (see Fig. [3|) that corresponds to an 
infinitesimal displacement field. As described by Eq. (jlip . tensor F describes the transformation 
from the relaxed, reference state (sphere, generic point X) to the current state (ellipsoid, generic 
point x). F is chosen to be a pure deformation, while dt may include some rotation. In 
this example, the evolution of tensor F is purely kinematic (see Appendix |A|) . as the material is 
supposed to behave in a reversible, elastic manner here. See Fig. [29] for details on the relaxed state. 

Let us now construct a consistent framework for the material elasticity, with the help of 
Figure [TOl which pictures the evolution of the stored deformation due to an applied flow. 

A. Stored deformation 

In such a plastic material, as pointed out above, there is no point to define the deformation 
with respect to some remote, initial reference state. At any time, however, every region of 
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the material is stretched (or not). In other words, its conformation differs (or not) from 
the conformation it would have in the absence of stress from the neighbouring regions: the 
stress locally induces {via its elasticity) a stored deformation, which can be visualized in the 
form of an ellipse (or, more generally, of an ellipsoid). To obtain the ellipse (see Figure [3D, 
one needs to cut out a piece of material, draw a circle on it while it is thus relaxed, and put 
it back in place. 

The circle (or, more generally, the sphere) in the relaxed state is described by: 



(10) 



where X is a vector whose origin is the center of the sphere, and whose end is a generic 
point on the sphere (see Figure [TOl) . The ellipsoid in the stretched material can be described 
by some vector x, which depends linearly on X if the sphere radius R is infinitesimal, and 
can be expressed in terms of the transformation from the relaxed state to the current state: 



x = F-X 



fill 



The transformation from X to x (see Figure [TOl) can be chosen as a pure deformation (see 
Fig. E]). Tensor F is thus symmetric, and in its principal axes, we have: 



^ Xi 0^ 



A2 
A3 



where, in the case of an incompressible material, detF = A1A2A3 = 1. 
The equation for the ellipsoid is obtained from equations ffTOl) and ffTTl) : 



(12) 



x^ ■ F'^ -x 



R' 



(13) 



B. Finger tensor and associated deformation 

As mentioned in the caption of Figure [3l the relaxed state local orientation can be chosen 
arbitrarily. This choice thus must not have any incidence on the material elasticity. We 
therefore need a variable that reflects local deformation without being sensitive to local 
orientation. A commonly used such variable is the Eulerian Finger tensor B = F ■ F^ . In 
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the present case where the local orientation is chosen such that F be a pure deformation 
{i.e., represented by a symmetric tensor), 



^ f3i 0^ 



B = F' 



p2 

/?3 



(14) 



where (3i = Xj. From the Finger tensor, we can also construct a deformation 48|| e = ^{B — T): 





|(/32 - 1) 

i(/33-i)y 



(15) 



C. Evolution of the stored deformation 

Let us now describe how a piece of material under stress is further deformed when an 
infinitesimal displacement field vdt is applied to the material. The corresponding deforma- 
tion, 

dv ■ 

Vw dt, where (V'u)ij = (16) 



weakly deforms the ellipsoid (see Figure [TUl) . As shown in Appendix |Xl one can express the 
equation for the deformed ellipsoid and derive the evolution equation for tensor B: 

dB 



dt 



-Vv ■ B - B - Vv 







(17) 



where d/dt = + ■u ■ V is the particulate derivative and where the entire left-hand side 
is the upper-convected objective derivative of B. 

The evolution of the associated deformation, defined by Eq. (ITSl) . also involves its upper- 
convected objective derivative: 



de 

wv ■ e — e ■ Vv = D, 

dt 



(18) 



where D is the applied deformation rate (see Eq. [2]). 

Equations f[T7l) and (ITS]) provide the variation of the stored deformation in the case of 
a purely elastic behaviour, i.e., in the absence of any plasticity. We shall soon discuss the 
origin of plasticity and indicate how it may enter such evolution equations. But let us first 
give a clear formulation of how the stress relates to the stored deformation. 
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D. Elasticity 



Since the material is supposed isotropic, the principal axes of the stress coincide with 



those of the Finger tensor, and the stress can be expressed as a linear combination 
three powers of B, for instance: 

a = aol + ai B + a2 B'^ 



49j of 



(19) 



where oq, ai and 02 are scalar, isotropic functions of B {i.e., of its scalar invariants). Equiv- 
alently, the deviatoric (traceless) part of the stress 



a = dev{a) = a — ^^^l*^) = o' + pi 
(where p is the pressure and d is the dimension of space) can be expressed as 



(20) 










o\ 


a = ai dev(5) + a2 dev(5^) = 

















^3) 



(21) 



and the pressure is given by: 



P=-ao~ ^tT{B) - ^tT{B') 



(22) 



More specifically, we will assume that the material is hyperelastic, i.e., that the stress 
results from the differentiation of an elastic potential E{B), here defined as being an elastic 
energy per unit mass of the material. It can then be expressed 3J] as: 



with: 



dE 

a = 2pdev -.5 



dE^ao , a, 02 
^ d5 2 2 2 



(23) 
(24) 

(25) 



In the case of an incompressible material (see Section IIII Gl for examples), the stored 
deformation (expressed in terms of B or e) does not depend on pressure: it only depends 
on the deviatoric part a of the stress. Conversely, the deviatoric stress a can be expressed 
in terms of -B or e: 

a = J^{B) = g{e) (26) 
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In this incompressible limit, the pressure p (more precisely, its term ao) varies very 
strongly when det B departs from unity. In practice, one usually considers that p is not 
known explicitely: in a practical situation, it must be obtained from the boundary condi- 
tions, while strictly enforcing the constraint det i? = 1. 



VI. PLASTICITY 



Now that we have described the evolution of the stored deformation in the material 
(Eq. flTTI) or flTSl) ) and the material elasticity (Section [V DP . let us turn to the description of 
its plastic properties. 

As mentioned earlier, a plastic term must be added to Eq. (ITTl) or (ITSl) in order to reflect 
the deformation due to topological events: 

'}^-Vv-B-B-Vv^ = -2D^ (27) 
dt P ^ ^ 

de 

— - Vv-e-e-Vv^ = D - (28) 

dt P ^ ^ 

This equation is the tensorial version of the decomposition of the applied deformation into 
a kinematic part and a creep part expressed by Eq. (jS]). 

In the present section, we recall the microscopic origin and discuss the mathematical 
properties of the plastic term in the evolution equation. 



A. Threshold for a single Tl process in two dimensions 

Figure [TT] provides a two-dimensional illustration of a Tl process in an emulsion or a foam. 
Four objects (hereafter called quadruplet) are initially organized in a diamond conformation. 
As compared to the initial conformation at rest (point Ri), a weak applied stress induces an 
elastic deformation, symbolized by a continuous branch (thick line) . Once the applied stress 
exceeds a certain value (symbolized by ay on Fig. [TTl) . the system undergoes a discrete flip 
while the four objects swap neighbours. This is the Tl process, symbolized by a jump onto 
the other branch. 

More precisely, let us suppose that the quadruplet is symmetrical like that on Fig. [TT], 
with axes ^ and (. The stress orientation that is most favourable for the Tl process to occur 
is that in which the stress axes are aligned with those of the quadruplet. Thus, Fig. [TT] is 
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FIG. 11: Schematized Tl process in two dimensions: conformations of four neighbouring objects 
(quadruplet) subjected to external forces. One pair of objects is oriented along axis ^ and the 
other one along axis ( (see inset). The stress principal axes are assumed to be aligned with those 
of the quadruplet, which therefore adopts a diamond configuration. Moreover, the applied pressure 
is assumed to be constant. The diagramme indicates how both pairs compare as for the (centre- 
to-centre) inter-object distance (horizontal axis) and as for the external stress that they undergo 
(vertical axis). By convention, positive values of the stress (resp. S^^) denote traction on the 
pair that is oriented along axis ^ (resp. Q. A Tl process starts from the relaxed situation Ri, 
denoted by a black circle. Through compression along axis ^ {i.e., < 0) or traction along ^ 
(S^^ > 0) or through a combination of both, the system follows branch 1 (thicker line) towards 
the critical point Ci. The Tl itself is the sudden transition from Ci (on branch 1) to branch 2 
(thinner line). In the case of an incompressible material the quadruplet conformation depends only 
on — S^^, while for a compressible material it additionally depends on -|- = —p. 
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drawn in terms of E^^ and S^^, with the non-diagonal stress component E^^ being equal to 
zero. 

The condition for the quadruplet to remain on branch 1 of Figure [TT] can be expressed in 
the form: 

/eigcn(E) < (29) 

where /eigen(E) = S^^ - - ay^^ (30) 

Here, ay is a real, positive number. When the pressure is varied or when the stress is not 
aligned with the quadruplet (non-zero shear stress component, i.e., 7^ in the present, 
2d case), it more generally takes the form: 

a/i = a/i(p,E^c) (31) 

ay^^ is expected to increase with the applied pressure p = — E^^ "^CC- Indeed, increasing 
the pressure reduces the typical object size, thus enhancing surface tension effects. It should 
therefore increase the stress threshold in such surface tension sensitive systems as foams or 
emulsions. 

In the case of an incompressible material, cr^^^ depends only on E^^, not on p. More 
generally (in two and three dimensions), the stress threshold function /eigen depends only on 
the deviatoric part of the stress and can be written in the form: 

/eigen (E) = fifcigen [deV ( E ) ] (32) 

As schematized on Fig. [T21 let us now discuss how the Tl-threshold discussed above 
impacts the plasticity threshold and the plastic flow. 

B. Plasticity threshold 

Instead of a single Tl-capable quadruplet (Fig. [121 left), let us now consider a mesoscopic 
element of material, containing quadruplets with many different orientations (Fig. [121 cen- 
tre). 

If we neglect the effects of disorder (see next paragraph for a short discussion), the plastic 
threshold for this chunk of material is reached when the most favourably oriented quadruplet 
reaches its own threshold: 

f{a) = max /eigen (0"a) 

a. 
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FIG. 12: From Tl process to plastic flow. Left: the stress needed for four bubbles to undergo 
a Tl process is illustrated on Fig. [TTJ Centre: when a chunk of material contains a number of 
quadruplets of bubbles with various orientations, it may be postulated that the plasticity threshold 
corresponds to the most favourably oriented quadruplet (this assumption discards the effects of 
non-homogeneous deformations, and is discussed in the main text). Right: the rate at which Tl 
processes occur in the chunk of material when it is subjected to elongation B is expressed in the 
form of a plastic flow , where the material configuration at time t (dashed line) is mapped onto 
configuration at time t + dt (full line) through infinitesimal deformation D^dt. 

where a runs over all orientations in space, and where cTq, is the representation of tensor a in 
the axes of a quadruplet oriented according to a. Note that the most sensitive quadruplets 
are those oriented along the stress: the most favourable orientation a coincides with that of 
the stress tensor (possibly up to some permutation of the axes). 

In the case of the particular, two-dimensional threshold given by Eqns. (l3T!l and (!30|) . we 
obtain: 

fi(^) = W{1)- (^{2)1- (Tyip) (33) 

where cr(i) and (T(2) are the eigenvalues of the stress tensor and where p = — o"(2) is the 
pressure. 

In the case of an incompressible material, the plastic threshold depends on the sole 
deviatoric part of the stress, like in Eq. (!32l) : 

/(cr) = maxfifeigen[dev(cra)] 

a 

= 9[devia)] (34) 
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C. Disorder and plasticity threshold 



The local disorder of soft object positions and interactions implies that the stress acting on 
a particular quadruplet slightly differs from the ambient mesoscopic stress, and the threshold 
of some of the quadruplets is lower than expected for a particular orientation of the applied 
stress. 

As a result, the threshold value Cy in Eq. fl33|) above is in fact slightly lower than the 
microscopic value 0"^^^ in Eq. (!30|) . Similarly, expression (!33|) for /(c) is slightly over- 
estimated: (Jq, should in fact be understood as the a-oriented representation of a locally 
disorder-enhanced version of the mesoscopic stress cr. 



D. Plastic deformation rate 



In a continuum model, the rate at which topological events occur is expressed as the 
plastic flow and depends on the stress a, or equivalently on the material deformation 
(tensor B) via its elasticity (see Eq. [2^ . The plastic flow therefore has the following form: 



hiB), 



(35) 



where h is a tensorial function. Since the plastic flow results from Tl events, whose principal 
axes almost 50|] coincide with those of the stress (see paragraphs IVl Al and IVl Bl) , function 



h is an isotropic function of B: it is such that the principal axes of also coincide with 
those of B. Thus, in the same axes as those of B: 

fx n n \ 



Di 



5i 

^2 
^3 



(36) 



Equivalently, since powers of tensor B constitute a basis for symmetric tensors that have 



the same principal axes, can be decomposed for instance in the following way: 



= bol + hB + b2 B\ 



(37) 



where scalar functions bo{B), bi{B) and 62 (-B) depend solely on the scalar invariants of 
tensor B. Let us emphasize the fact that in general, all three functions should be non-zero 
above the plasticity threshold. Restricting to be proportional to B, for instance, would 
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mean that all three directions would flow indepedently of each other. This would be a very 
specific choice, certainly not relevant for most materials. It would be the plastic equivalent 
to taking a Poisson ratio u = in linear elasticity, a property valid for only a restricted class 
of materials. 

Function h (or equivalently eigenvalues Si in Eq. [Snior coefficients Bq, bi and 62 in Eq. [HTIl . 
must also obey some constraints based on physical grounds, which are discussed in the next 
few paragraphs. 

E. Incompressible plastic flow rate 




D dt = dB + dt 
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FIG. 13: Incompressibility and plasticity. The incompressible evolution Dp of the relaxed confor- 
mation implies the usual property tr(Z)^) = 0. Despite this, the corresponding flow in real 
space is not traceless in general (tr(L'^) / 0): it is related to through a finite transformation 
B, see Eq. p8]) for details. 



Due to plasticity, the deformation between two real configurations (see Figure [12]) does 
not refiect solely the increment ^dt in stored elastic deformation, but also the irreversible 
drift Dpdt of the relaxed local configuration, where the drift rate is related to the plastic 
term defined by Eq. fl35l) through covariant transport via tensor F: 



(38) 
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In the present work, we do not consider all possible plastic deformation modes. Among 
the situations depicted on Fig. [HI conformations (0), (1) and (4) are equilibrium situations 
(at least on time scales where the dispersed phase does not diffuse and where the continuous 
phase does not permeate). Hence, modes (0) (1) and (0) (4) correspond to plastic 
evolutions of the material. In the present work, as stated in Section IIII F 11 such evolutions 
that involve mutual diffusion between both phases are not addressed. 

The plastic evolutions of interest for us correspond to stress-induced evolutions, where 
the relaxed conformation is always characterized by the same amount of material per bubble 
or droplet (inside and around it). Among the equilibrium conformations (0), (1) and (4), 
only conformation (0) is thus eligible. 

We are thus interested in elastic deformations (combinations of modes (0) (2) and 
(0) (3)) followed by plastic rearrangements that bring the local conformation back to 
(or at least towards) situation (0). In other words, the local relaxed conformation is always 
locally similar to situation (0), even though some bubbles or droplets may have swapped 
positions. 

In particular, the evolution of the local relaxed state must not be accompanied by any 
change in volume. This can be expressed through a condition on the drift rate D^: 

tr{D^) = (39) 

Using Eq. (|38|1 . this can be expressed in terms of D^: 

= tr(5-^ ■ D^) (40) 

- _L I ^ 

" /?1 /32 /^S 

Expression (1401) thus expresses incompressibility for the plastic flow rate . 

Note that in order to satisfy the incompressibility condition given by Eq. (l40i) . the plastic 
flow given by Eq. ( !37j) can be written in the form: 

= biB -deviB) + b2B -deviB^) (41) 

Scalar functions bi and 62, hke 60, ^1 and 62 in Eq- (1371) . depend solely on the invariants of 
tensor B. 
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In this case, the coefficients of Eq. (1371) can be obtained through the relations: 



62 = &i + b2 tr(5) 



h = -|[tr(5)]^ + |tr(5^)-|tr(5) 



bo = b2 



(42) 
(43) 
(44) 



and are related to one another through: 



([tr(5)]2-tr(i?2))+6, + |tr(i?)=0 



(45) 



F. Elastic versus plastic incompressibility 

At this point, it may be useful to precisely delineate two different types of incompress- 
ibility, which can be discussed using Fig. [HI 

As stated in paragraph IVIEI above, in the present work we are only interested in locally 
relaxed states of type (0), and the plastic flows we consider correspond to relaxed conforma- 
tions that evolve among such states of type (0), and the density of the locally relaxed state 
is conserved: the material is plastically incompressible. 

By contrast, as mentioned in paragraph IIIIF 11 the actual state of the system (with 
local stresses) may have a local density that differs from that of the relaxed state, due to 
evolutions of type (0) (2). In other words, the material is assumed to be compressible, 
that is: elastically compressible.\51\ 

It follows that with our assumptions, the material may locally change volume (and the 
density p then departs from its initial value po) when the flow is not divergence-free (V • 
V = trD 7^ 0). Nevertheless, because the plastic evolution is assumed incompressible, the 
material returns to its initial density po as soon as the local stress vanishes. 

For some systems, elastic compressibility is negligible and the material can be considered 
elastically incompressible, as explained in paragraph IIII Gl The material then locally never 
changes volume, and the flow is then divergence- free (V ■ v = tiD = 0). 

G. Thermodynamic constraints 

The total work developed by the stress in the material is given by tr(o" ■ D). In the 
present system, the internal energy U is purely elastic: U = E. The first principle of 
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thermodynamics can thus be written as 

pdE = tr(a ■D)dt + 6Q, (46) 

where 6Q is the heat uptake per unit volume. 

Besides, the dissipated power Pdissip per unit volume results from viscosity, especially 
during the relaxation of individual Tl processes. It constitutes the only source of entropy 
in the system, which is non-negative according to the second principle: 



dissip 



TdS-crcatcd >0 (47) 



It also constitutes the only source of heat in the system. If we assume that the system 
remains at a constant temperature, this condition can be written as: 

-Pdissip 

dt + 6Q = (48) 

The Clausius-Duhem inequality is readily derived from the above equations: 

dP , , 

tr(cT ■D)-p — >0 (49) 

Let us express this inequality in a different way. 
From Eq. (l23l) . we get: 

tr((7 ■ D) = 2pti (^B (50) 

We also have: 

dE fdE dB\ 

where dB/dt is given by Eq. fl27|) . Hence, from Eqs. fj49|) . fl50|) and flSTl) . we obtain the 
constraint on the plastic flow that corresponds to the Clausius-Duhem inequality: 

Pdissip = 2p tr (d^ ■ ^) > (52) 



H. Plastic deformation rate: a summary 

In the last few paragraphs, we showed that the plastic deformation rate must obey con- 
straints that express the fact that: 
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1. the plastic flow depends on the stress or equivalently on the stored deformation (see 
Eqs.[35]or[37D; 

2. if applicable, the plastic flow is incompressible (see Eqs. HOlor I^Tjl: 

3. the associated dissipation is positive (see Eq. 1321) . 

VII. COMPLETE, CONTINUOUS MODEL 

In this brief section, let us discuss how the constitutive equation derived in the present 
work can be inserted into a set of equations and provide a complete, continuous model. 

Keeping in mind the considerations of paragraph IIIIF[ let us now close the evolution 
equation of the system (Eqs. [27| or 1251) 

d-B 

— -Vv- B - B -VtF = -2D^ (53) 
dt P ^ ^ 

One needs the elastic law (Eq. [19] or [23]) 

dE , , 

a = 2p—-B (54) 

and the usual force balance equation 

-> dv 

V-a + p/ = p-, (55) 

where / represents external forces per unit mass. The evolution of the density p obeys the 
usual mass conservation equation 

|^ + V-(pt;) = ^ + ptrD = (56) 

We provide later (see Eq. [SI]) a version of this equation that includes liquid permeation. 

Let us now mention the special case of an initially non-dilated material, and discuss 
dissipation at weak applied stresses. 

A. Initially non-dilated material 

Let us assume that the material initially has a uniform density i.e., with stored deforma- 
tion 

p{to,r} = po (57) 
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and that it verifies everywhere 

detS(to,r) = 1, (58) 

In such a case, our assumptions stated in paragraph llll F ll imply that the material density 
is related to the determinant of tensor i? in a very simply way. Indeed, the derivation shown 
in Appendix [B] implies that at any later time, throughout the sample: 

P = (59) 

With such assumptions, Eq. (15^ can therefore be used to replace p within Eq. (15^ . and 
Eq. (!56|) is then not useful any more. 

B. Viscous losses under weak stresses 

The Bingham model addressed in the present work reduces to a simple elastic system when 
subjected to weak stresses (below the yield stress). As a result, vibrations may be present in 
the material, at arbitrary high frequencies. Such vibrations may be undesirable. Not only do 
they make numerical simulations of the above system of equations more complicated, if not 
impossible, but they do not faithfully reflect the damping observed in real materials. This 
problem does not arise under large stresses, as the plastic flow rate introduces dissipation. 

One of the simplest ways to introduce some dissipation at weak stresses is to add a viscous 
term to the stress. Eq. (l54l) thus becomes: 

(T = 2p^-B + V{D, B) (60) 



where the viscous term is a symmetric tensor which depends linearly[52|] on the applied 
deformation rate D. As the material is locally anisotropic due to the stored deformation, V 
additionally depends on B. 

The general form of function V{D, B) is a sum of terms whose principal axes are those 
of tensor B: 

h{B)tT{l2{B)-D) (61) 
and of terms that depend tensorially on D: 

mi{B) ■ D ■ m2{B) + m2{B) ■ D ■ mi{B) (62) 

Here, all functions Zj and rrii are isotropic scalar functions of tensor B. 
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If we neglect any effect of the material deformation, we may simply add a linear viscous 
term, and Eq. fl5^ then becomes: 

dE 

a = 2p—-B + riD (63) 

VIII. HOMOGENEOUS CONSTANT SHEAR FLOW: ONE EXAMPLE 

Let us now use the constitutive Eq. fl27j) to obtain the evolution of the system in a 
very common type of flow: a shear flow at a constant shear rate 7 starting at time t = 0. 
For non- homogeneous flows, the constitutive equation must be combined with the classical 
mechanics equations for continuum media, as mentioned in Section IVlIi Here, for simplicity, 
we assume that the material is homogeneous and incompressible (uniform density po all 
times) and that the flow remains homogeneous: no shear-banding, etc. 

A typical result from Eq. (!27|) . coupled with an elastic law (here incompressible, see 
Eq. [2T]) is the shear stress as a function of the shear deformation 7 since time t = when shear 
started. An example of such a mechanical response (with parameters chosen as described in 
Section rVIIIBI) is provided on Fig. [141 The equations of our model provide: 

1. the material response in its elastic state (first part of the curves on Fig. [Ml) : 

2. the threshold that marks the onset of plasticity (point where curves split apart de- 
pending on the shear rate); 

3. the transient response that results from plasticity (second part of the curves); 

4. the stationary response, as a function of the applied shear rate 7 (see Fig. [T6l) . 



A. Method 

Let us recall the evolution equation ([27]) for the material deformation (tensor B): 

d /? 

— -Vv- B - B -Vif = -2D^ (64) 
dt P ^ ^ 

This equation must be understood as written in a basis attached to the (fixed) laboratory 
frame. Let x be the axis of velocity, y the axis of the velocity gradient and z the vorticity 
axis. The symmetry of the shear flow implies that the corresponding material deformation 
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FIG. 14: Shear stress ai2 in the course of a shear experiment, as a function of the shear deformation 
7, for three different shear rate values 7. The elasticity and plasticity terms are chosen as described 
in Section IVIIIBI The response increases during the initial period of time, when the material 
deforms purely elastically (as expected, the response is then independent of the shear rate). The 
increase then levels off once the plasticity threshold has been reached, and saturates at a stationary 
value. For a plot of the final, stationary value on the applied shear rate, see Fig. [161 

B has two principal axes within the xy plane (in directions X and Y yet to be determined, 
see Fig. [T7j) and one principal axis along z. Note that and a also has the same principal 
axes X, Y and z as B (see Eqs. [251 and I^Tl) . 

As shown in Appendix I C 1\ Eq. fl64p provides a system of differential equations for 
and for angle 6 between and y (see Fig. [TT]) : 

u = cos(2^) (65) 

^ = 7/5iv^r^ - 26i{Pu P2) (66) 
d/?2 



-7/92 - 262{(3u P2) (67) 

dw . n K { . /3l + l32 



where Eq. fl69|l results from the assumed material incompresibility. 
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FIG. 15: Dissipation Pdissip psr unit volume and injected power ai2 7 in the course of a shear 
experiment, as a function of the shear deformation 7. The apphed shear rate is 7 = 8. The 
elasticity and plasticity terms are chosen as described in Section IVIIIBI Note that there is no 
dissipation (Pdissip = 0) in the elastic regime, prior to the onset of plasticity. Part of the injected 
work is not dissipated (non-zero surface area between both curves) and corresponds to the elastic 
energy stored in the material. 



The above equations provide the evolution of the material deformation (tensor B) from 



the initial situation at rest 



53| to the onset of plasticity, the plastic transient and the final, 



stationary state (see Section [(T4l of Appendix). 

Once the evolution of tensor B is known, the stress is obtained through Eq. ([2] 
As for dissipation, given by Eq. (!52|) . it involves both the gradient of the elastic energy 

(see Eq. [25]) and the plastic deformation rate given by Eq. fj^T]) . 
The dissipation per unit volume can now be expressed as: 



dissip 



dB 



(70) 



+ (02^0 + ctibi) ti{B) 
+ {aib2 + a2hi) ti{B^) 
+0262 tr(5=^) 



(71) 
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FIG. 16: Shear stress crl2 and dissipation power -Pdissip p6r unit volume in the stationary regime, 
as a function of the shear rate 7. The elasticity and plasticity terms are chosen as described in 
Section Emll 



y 




X 



FIG. 17: The homogeneous shear flow imposed in basis xyz attached to the laboratory [z being 
the direction of vorticity) causes the principal axes of tensor B to tilt (basis XY z) with respect to 
the shear. The lengths of the ellipse axes represent the magnitudes of eigenvalues j3i and [32- 
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When Dp is known explicitely, it can be calculated more directly: 

-Pdissip = (ai + a2/3i) 6i 

+ (fli + 02/92) ^2 

+ (ai + 02/93)53 (72) 

Figure [15] displays both the injected power, al2 7, and the dissipated power, Pdissip, in 
the course of a shear experiment. Notice that as long as the material remains in the elastic 
regime, no dissipation occurs. The injected power is being stored entirely as elastic energy. 
As dissipation starts, an overshoot of dissipated power is observed. Asymptotically, both 
quantities converge towards the same value while the elastic deformation and energy of the 
material reach their stationary values. 



B. Chosen elasticity, threshold and plastic flow 



Let us choose a very common (Mooney-Rivlin) type of incompressible elasticity, which 
has been shown to suitably approximate the non-linear elasticity of foams 35|, l36|]. The 
corresponding elastic energy |3J| can be expressed as: 



Poi? = |(lB-3) + |(ll^-3) 



(73) 



where 



Ib = ti{B) 

11b = ^[ti'iB)-ti{B')]=tiiB-') 



Correspondingly, the coefficients of Eq. (I2T]) can be expressed as: 



(74) 
(75) 



fli = ki + k2lB 

0-2 = -k2 



As in refs. 



35 



(76) 
(77) 



36(1, we choose the values of ki and k2 in terms of the shear modulus G as: 



k. = \G 



(78) 
(79) 
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Note that in Section IIXI we explore other values for ki and ^2 (still in the framework of a 
Mooney-Rivlin elasticity). 

For the plastic deformation rate , as expressed by Eq. (HTi) so as to obey plastic 
incompressibility, we choose the following coefficient values: 

bi = [iiB - 4) eitiB - 4) (80) 
62 = (81) 

where 9{x) = 1 when x > and 9{x) = otherwise. Note that in Section HXl we explore 
the case 62 = ^1 for comparison. 

This choice implies, in particular, that the plasticity threshold corresponds to the follow- 
ing condition: 

tr(5 - 4) = (82) 

We now discuss qualitatively the evolution of tensor B. We then show the effects of the 
material elasticity and discuss the corresponding rheological response. 



C. Three-dimensional evolution of the material stored deformation 

With the choices made in Paragraph IVIIIBI for the material elasticity and plasticity, the 
evolution of the system given by Eqs. (I65|) - (l69|) is depicted on Fig. (fTSl) . For three different 
values of the shear rate 7, it represents the trajectory of the material stored deformation as 
for (3i and (32, successively in the elastic regime and in the plastic regime. The plasticity 
threshold, given by Eq. fl82|) is also represented, as well as the locus of the stationary states. 

In the elastic regime, the form of Eqs. (lC2p - (lC3p imply that /^a remains equal to its initial 
value, = 1, i.e., that the product P1P2 remains equal to unity. This can indeed be checked 
from Eqs. flCll|) - flC12|) . and it reflects the fact that as long as no plastic events have occured 



in the material, the absence of material deformation in direction z (planar shear in the x-y 
plane, see Eq. ICip implies that the stored deformation is not modified in direction z. 

In other words, in the elastic regime, the fact that no plasticity occurs in direction z {i.e., 
63 = 0) implies that in the same direction, the stored deformation remains constant, i.e., 
P3 = const. 

Conversely, in the stationary regime, /S^ remains constant (like most other quantities) 
and this imposes that no plasticity occurs in direction z, i.e., 63 = 0. 
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Meanwhile, in the transient regime, the stored deformation (3^ in direction z has evolved 
from its initial value towards its new, stationary value, despite the absence of any velocity 
gradient in this direction. This is made possible by the plastic events (^3 7^ 0) which allow 
internal relaxation within the material even in the absence of any flow in this direction. 
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FIG. 18: Form of the stored elastic deformation in the com'se of an experiment and in the stationary 
regime, for three different values of the shear rate 7. The axes are the first two eigenvalues, /3i and 
/32, of tensor B. 



D. Shear thinning 



The stationary shear stress is represented on Fig. [11] as a function of the shear rate 7. 
Above the yield stress, the stress increases when 7 is increased, as expected. Notice that 
it increases in a subliminear way. We believe that the main reason for this shear-thinning 
behaviour could be the fact that the plastic flow rate we chose, as expressed in termes 5, is 
cubic rather than linear for large stored deformations: 



= (trB - 4) 9{tTB - 4) 



- — tr(fi) 



(83) 
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FIG. 19: Eigenvalues of the stored elastic deformation in the course of an experiment as a function 
of the shear deformation 7. The applied shear rate is 7 = 8. 
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FIG. 20: Orientation of the stored elastic deformation in the course of an experiment as a function 
of the shear deformation 7. The applied shear rate is 7 = 8. 

where, again, 9{x) = 1 when a; > and 9{x) = otherwise. 
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E. Normal stress differences 



The stress tensor a is obtained from the Finger tensor B through Eq. (!2T]) . In the case 
of a homogeneous shear flow (see above), it can be expressed not only in the XYz basis 
(Eq. [?T]) but also in the xyz basis associated with the flow: 



a 



0'12 cr22 
(73 



cs{ai — s^ai + c^a2 
as 



(84) 



where s (resp., c) denotes the sine (resp., the cosine) of angle 9 between axes x and X 
(see Fig.[I7D- 

The first and second normal stress difference can be expressed as: 



Ni = au - 0-22 
= (o-i - 0-2)^ 

-^2 = 0-22 — 0-3 

1-u _ 1+u _ 
= - 0-3) H ^^(^^2 - (^3 J 



(85) 
(86) 



with (Tj given by Eq. fl2T]) : 



CTi — CT'i 



ai (A - Pj) + a2 {(3f - P') 



(87) 



Hence, for instance: 



iVi = u[a^{(3,-(32) + a2{f3f-(3i)] 



No 



+ 



/5i + /32--^ + w(/32-/3i) 
P1P2 



^2 



(A/52) 



+ w(/?2-/3?) 



(88) 



(89) 



In the case of Mooney-Rivlin, coefficients ai and a2 are given by Eq. (1761) and the normal 
stress differences can be expressed as: 



47 





FIG. 22: Shear stress ai2, deviatoric stress in the vorticity direction and normal stress differences 
Ni and in the stationary regime, as a function of the shear rate 7. 
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N2 



- (32) 



h + 



P1P2 



(90) 



2 

k 

+ - 



P1 + P2- 



+ U{^2 - A) 



2piP2 V u- 



Pl(32 



P1P2 



(91) 



F. Discussion of the stress response 



Let us now comment briefly on the results presented on Fig. 
Normal stress differences A'^i and increase gently and monotonically to reach their 
stationary values. 

The material is under traction in the direction of vorticity ((T3 > 0), as well as in the 
direction of while it is under compression in the direction of ^2 (not shown). This is 
consistent with the fact that it is stretched in the direction of vorticity {/J^ > 1). 

The salient feature of these results is that the shear stress presents an overshoot 
during its transient, a behaviour which is intrinsically unstable in a homogeneous material, 
and could trigger flow localization. Such an overshoot for the shear stress has been observed 
in sheared foams 



IX. IMPACT OF ELASTICITY AND PLASTICITY ON THE STRESS RE- 
SPONSE 

Our dynamical equation was formulated in terms of tensor B (see Eq. [MI) . From the 
resulting evolution of B, when can then derive the stress evolution using the elastic law 
(Eq. [19]). In this section, we explore how the choice for the elasticity and the plasticity 
impact the stress response. 



A. Impact of elasticity 

Note that in this paragraph the deformation history will be the same in all cases, and 
only the stress in the material will differ. Just like in Section IVIII| we restrict ourselves to 
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a Mooney-Rivlin type of elasticity, which is complex enough for illustrative purposes, but 
which has no specific properties with regards to our problem. 



The elasticity of dry foams has been shown 
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36| to be well captured by such a Mooney- 



Rivlin elasticity, using parameter values fci = 1/7 and k2 = 6/7. These values correspond 
to Figs. [T1H22] where we took the modulus value G = 1. 

Exploring different values for ki and k2 (still with G = 1), we show in Fig. [23] the time 
evolution of shear stress, stress in the vorticity direction and normal stress differences. The 
mechanical behaviour remains similar in all these cases. In particular, the shear stress always 
presents an overshoot. This overshoot is more pronounced for ki = 1 and k2 = 0. In all cases, 
the normal stress differences do not present any overshoot and reach monotonically their 
stationary values. As for the stress along the vorticity direction. It presents an undershoot 
when ki > k2, and its stationary value becomes negative at some point between the situation 
in the center of Fig. [23] {ki = k2 = 1/2) and the situation on the right-hand side {ki = 1 
and k2 = 0). 



B. Impact of plasticity 

Let us now briefly explore the impact of the plastic deformation rate D^{B) through the 
influence of parameter 62 in Eq. ( 1411) . 

In Fig. [Ml we take hi = {trB - 4)6'(trfi - 4) like in Figs. [HHZS We then observe the 
(unspectacular) effect of choosing 62 = ^1 rather than 62 = 0. Note that in this case, as 
opposed to the role of elasticity in paragraph IIX A| the time evolution of tensor B is affected 
by the choice of the plastic law. 



C. Discussion 



In paragraphs IIX Al and IIXBI above, we illustrated the fact that the stress response of 
the material can be altered through a change in the material elasticity and plasticity alone. 

Our approach emphasizes the fact that the choice of the convective derivative is some- 
what arbitrary. A material is often considered as rather "upper-convected" or rather "lower- 
convected" because its stress evolution follows more closely the eponymous convective deriva- 
tive. In fact, always choosing the upper-convected derivative (and tensor B for a measure 
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of the deformation) would be equivalent, provided the elastic and plastic laws be suitably 
adjusted. A forthcoming paper 32| will discuss in detail the consistency of commonly used 
rheological laws (see Section IIVCI above) and explore these issues more in depth. 



X. CONCLUSION AND PERSPECTIVES 

In order to describe the rheological behaviour of foams and emulsions (and potentially 
other materials), we have developped a continuous framework to describe the evolution of 
the stored deformation tensor. It includes elasticity, up to the large deformations commonly 
encountered in such systems, and plasticity. 

Within this framework, we showed that it is possible to play with various expressions for 
the elasticity and the plasticity. We hope that it is thus possible to adequately describe the 
rheological behaviour of a large range of foams and emulsions with dispersed phase volume 
fraction close to unity. 

In less concentrated foams and emulsions, volume fraction should be added as an extra 
field to account for osmotic compressibility of the dispersed phase. The issue whether the 
model may be adapted to describe usual, aging complex fluids, is left for future investigations. 

As mentioned in paragraph IIII At the Burger model (see Fig. [25]) describes the response 
of dry foams under weak stresses adequately. By contrast, the Bingham model (used as 
a basis for our formulation, see Fig. [26ll captures the large stress behaviour. It would be 
interesting to combine both models, as represented in Fig. [271 We did not elaborate on this 
combined model, as it is technically more complex than our Bingham-like model, and we 
did not want to focus on short timescales or weak stresses, where both models have differing 
behaviours. 

A. Determining the elasticity from experiments 

Beyond all these rheological and phenomenological models, it is important to make the 
connection with the microscopic scale. 

From this point of view, the main issue is the definition of the deformation in the material. 
In the present work, deformation is built on a thought experiment: we cut out a piece of 
material at sufficiently large a length scale for disorder to be smoothed out. 
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FIG. 23: Impact of elasticity on the stress response. The shear stress CT12, the deviatoric stress 
in the vorticity direction and the normal stress differences A'^i and N2 are plotted as a function of 
the shear deformation 7, for a shear rate 7 = 8. Mooney-Rivlin elasticity is chosen with different 
parameter values. Left: ki = and ^2 = 1. Center: ki = 1/2 and k2 = 1/2. Right: ki = 1 and 
A;2 = 0. 
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FIG. 24: Impact of plasticity on the stress response. The shear stress cti2, the deviatoric stress 
in the vorticity direction and the normal stress differences A'^i and N2 are plotted as a function of 
the shear deformation 7, for a shear rate 7 = 8. Mooney-Rivlin elasticity is chosen with ki = and 
k2 = 1. The plastic deformation rate is taken as expressed by Eq. (j4ip with bi = (trS— 4)^(tri?— 4). 
Parameter 62 is varied. Left: 62 = 0. Right: 62 = 61. 

By contrast, a way to define deformation was introduced by Aubouy et al. |[25|] for all sys- 
tems in which the individual objects and their mutual contacts are experimentally accessible, 
such as 2D foams and emulsions (and 3D foams and emulsions when tomography will have 
become a routine technique for imaging such systems on Theologically relevant timescales). 



They first construct a symmetric tensor M from the centre-to-centre vectors [38|] for pairs of 
first neighbours. This tensor is a dilation (proportional to the unit tensor when the material 
is at rest). They then define the deformation U as the logarithm of M. 
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One might think that their tensor U is just one deformation among many others. But 
with respect to the particular systems that they consider (2D foams in the dry hmit under 
quasistatic deformation), it plays a very special role: it is the deformation for which the 



elastic law remains linear up to large deformations [26| (in fact, up to the onset of plasticity). 
This result confers some weight to the initial, scale-independence arguments 25|] for the 
choice of the logarithm as the link between their texture (dilation) tensor M and their 
deformation measure U = log M. 

In our perspective, as we will show in a more elaborate manner elsewhere 32|, the choice 
of a deformation measure corresponds to a choice for the convective derivative. The defor- 



mation measure U 



25 
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could therefore lead directly to a continuum formalism well-suited 



for (at least) two-dimensional foams in the elastic regime. 



B. Plasticity and mechanical noise 

In the present work, the plastic deformation rate was assumed to depend only on the 
local stress via the local stored deformation, see Eq. (!35l) . 

This dependence is sufficient to account for the non-local elastic effects observed in 



foams 



3, 



21| . and mentioned in paragraph II C[ Indeed, when plastic events occur (non-zero 



Dp), this impacts the local stored deformation via Eq. fl7r|) . This, in turn, affects the stress 
due to elasticity (Eq. 1231) . Equation (l55l) then implies that the stress is modified in the 
surrounding material. 

Thus, within this continuum model in which the plastic deformation rate depends solely 
on the stored deformation (or stress), plasticity at one location alters the stress (and the 
stored deformation) elsewhere in the material, thus possibly contributing to triggering plas- 
ticity there. 

Nevertheless, stress may not be the only factor that determines the rate at which Tl 
processes occur. For instance, in a recent work by Marmottant and Graner 2^], the plastic 
deformation rate {D^ in our notation) is proportional to the total deformation rate D when 
it has the same sign as the stored deformation 5J]. Incidentally, such a choice implies that 



relaxation 55|] cannot take place in the system. 



The fact that depends on D can be interpreted physically as the fact that mechanical 
noise may well help triggering plastic events. 
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FIG. 25: Burger model: schematic diagramme (left) and plot (right) of plastic flow as a function of 
tensor B, in the stationary regime and under small applied streses. In Burger's model, suitable for 
foams under low applied stresses (see paragraph IIII A|) . the plastic flow increases with the applied 
stress. If we neglect the short time scale response provided by the Voigt-Kelvin element (spring 
and viscous element in parallel), then this increase is weak and linear, as determined by the viscous 
element in series. 

As an alternative implementation of the impact of mechanical noise, one might take the 
plastic deformation rate as slightly enhanced in the presence of non-zero total defor- 
mation rate (although mainly determined by the stored deformation). For instance, if we 
assume that the mechanical noise is isotropic, we may include a multiplicative factor of the 
form (1 + r^) into , where plasticity would be significantly enhanced for deformation 
rates around and above r~^. 

We did not enter such subtelties in the present work, and restricted our study to a purely 
stress-dependent plastic deformation rate, as expressed by Eq. (I35l) . 

C. Plasticity and fluidity or fragility 

In the present model, the evolution of the stored deformation is given by Eq. (!27l) : 



Vv-B-B-Vv 

dt 




In the present work, the plastic flow rate 



was assumed to depend solely on the stored 
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FIG. 26: Bingham model: schematic diagramme (left) and plot (right) of plastic flow as a function 
of tensor B, in the stationary regime. In a model based on Bingham's model, such as that developed 
in the present work, the plastic flow is zero up to a certain stress (yield stress). Beyond the yield 
stress, it increases (linearly with the applied stress). 




B 



FIG. 27: Combined Bingham-Burger model: schematic diagramme (left) and plot (right) of plastic 
flow as a function of tensor B, in the stationary regime. In a combined Bingham-Burger model, the 
plastic flow would increase weakly with the applied stress at low stresses, and more strongly above 
the yield stress. We believe it would mimic the rheological behaviour of a foam quite adequately. 



56 



deformation B: 

= D^{B) (92) 

More generally, as mentioned in paragraph IX Bl above, one can include the effect of, say, 
mechanical noise by including a dependence of on the applied deformation rate D: 

D^ = D^iB,D) (93) 

In reality, the plastic flow rate may well depend on yet other variables than B and 

D. 

A few years ago, in their study of sheared two-dimensional foams between two solid plates, 
Kabla and Debregeas noticed that Tl events appeared preferentially in regions where the 



stress tensor was most disordered 



2l| 



Very recently, Goyon and coworkers 39|] studied the flow of three-dimensional emulsions 
in milli-fiuidic channels. Because of the large aspect ratio of the channel section, they 
were able to minimize edge effects and to obtain two-dimensional velocity profiles for the 
stationary flow of such emulsions between two walls. 

One of their key results is that even though they are able to measure locally the relation 
between the shear stress (which they derive from the applied pressure) and the shear rate 
(which they observe), they do not obtain a single mastercurve a = cr{'j) when changing the 
applied pressure. 

By contrast, Eq. ( 127|) . combined with Eq. ( 193|) . would predict, in stationary flow, that the 
stored deformation could be expressed in terms of the velocity gradient, B = B{S/v), and 
similarly for the stress, via elasticity. Hence, in the present situation, they would predict: 

a = (7(7) 



The observed behaviour 39|] contradicts this prediction. It thus shows that the plastic 
flow rate must depend on additional variables. The observations have been shown |39| to be 
compatible with the existence of diffusive, scalar quantity F, which they call fluidity. 

It will be a challenge, in future studies, to identify the microscopic origin of such a 
quantity, which should be truly tensorial in more general situations than plane shear. In the 
quasistatic limit, fluidity should probably be related to the frozen stress disorder identified 



by Kabla and Debregeas 



2l|. 
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D. Liquid permeation 




FIG. 28: Foam with inhomogeneities in liquid volume fraction. The gas/liquid interface is flat in 
inter-bubble films, and curved in vertex regions at the junction between 3 (in two dimensions) or 
4 (in three dimensions) bubbles. This curvature implies that the liquid pressure in the vertices is 
smaller than the gas pressure in the neighbouring bubbles. In wet regions, the interface curvature 
is less pronounced than in dry regions. As a result, there exists a pressure gradient, and the liquid 
tends to flow from wet regions towards dry regions. The intensity of the liquid flow depends on the 
hydraulic resistance in the Plateau borders that convey most of the liquid (depicted as the channel 
in this two-dimensional drawing). 

Permeation of the continuous, liquid phase through the network of channels formed by 
the bubbles or droplets may result from gravity, due to the density mismatch between both 
phases, and is called drainage. Even in the absence of gravity, permeation may occur when 
the volume fraction of both phases is not uniform throughout the foam. Indeed, as depicted 
on Fig. [2S1 the pressure in the liquid phase is lower in dryer regions than in wetter regions 
for otherwise identical pressure values in the gas bubbles. As a result, the liquid permeates 
from wetter regions towards dryer regions. 

Since vector field v represents the bubble (or droplet) velocity while the material density p 
includes the mass of the liquid phase (see paragraph IIII F\} . the density conservation Eq. (|56|) 
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must include a term that reflects permeation. 

This effect can be incorporated into Eq. fl56|) in a rough manner by adding a diffusion 
term: 

dp 

+ V ■{pv)= /^permeation Ap (94) 

where the diffusion coefficient -Dpermeation depends on the geometrical dimensions of the 
Plateau borders between foam vertices (which depend mainly on the volume fraction ip) but 
also on the hydrodynamic boundary conditions along the Plateau borders (which depend on 
various preparation conditions such as surfactant nature and concentration, salt, etc, in a 
very non-trivial manner). 



E. Towards dilatancy 

Eq. (IMI) above is only the first step towards a model that would include other effects 
known to exist in foams and emulsions. Indeed, relaxing the assumptions made in para- 
graph IIIIF II would enable us to include an important phenomenon which is well known in 
the context of granular media and which has been recently demonstrated in liquid foams 40] , 



namely dilatancy. When deformed, the local conformation of the foam, as schematized on 
Fig. [HI tends to go from conformation (0) to conformation (4). The physical origin of such a 
phenomenon is not clearly understood yet. It might be related to the thickening of films 41 1 



or Plateau borders 



41 



42l ] observed when a Plateau border glides at a solid wall. 
Together with the identification of the microscopic origin of fluidity (see paragraph IX CI 
above), including the effect of dilatancy will thus be yet another challenge in the ever- 
bewildering rheology of foams and emulsions. 
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FIG. 29: Augmented version of FigurefTOl As mentioned in the caption of Figure[3l the relaxed state 
local orientation is chosen in such a way that tensor F is a pure deformation. As a consequence, 
points X and X' do not coincide except for particular values of ViT. 

APPENDIX A: EVOLUTION OF FINGER TENSOR 

In order to derive the evolution of tensor B (see Eq. [T4l) in some known velocity field, we 
compute how an ellipse with generic point x (see Eq. [T3l) is convected into a new ellipse (see 
Fig. [TOl) . The generic point x' of the new ellipse is obtained as: 



il + Vvdt) -x 



(Al) 



The equation for the new ellipse can be written in two ways: 




Vv^ dt) . [F~\ . (1 



Vv dt) -x' = R^ 



The evolution of tensor B 



F-2 



is thus given by: 



d_ 
di 



[B-^] = -Vv^ ■ B 



-1 



B 



-1 



■ Vv 



(A2) 
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From Eq. (lA2p above, we then compute 



56] the evolution of tensor B, see Eq. (fT7|) . 



APPENDIX B: DENSITY AND STORED DEFORMATION 

Our assumptions stated in paragraph IIIIF II imply that the material density is related 
to the determinant of tensor i? in a very simply way. The evolution of det B can always be 
written in the form: 

d(det B) 
dt 

After multiplying Eq. fl27p by B^^, taking the trace and using the plastic incompressibility 
expressed by Eq. (1401) . we insert the result into Eq. (IBip and obtain: 

d(det B 



tr 



(deti.)i.-.f 



(Bl) 



dt 

and finally: 

d(l/v/ditS) 



2{detB) (tiD) (B2) 



, [1/VdetB) (tiD) = (B3) 

Comparing Eqs. (l56l) and (1B3P shows that for any given element of material, the current 
density and deformation are linked via their initial values: 



p{to, ro) ^ det B{t,r} 
where is the position of material point r at time to- 



p{t,r) _ ^/detB{to,ro) ^^^^ 



APPENDIX C: SHEAR FLOW 



1. Derivation of the shear flow equations 



In order to derive the evolution of tensor B (Eqs. [651468I) . let us use the notations of 
Fig. [13 

The shear velocity gradient, the material deformation B and the flow rate Dj^ have the 
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following form in basis xyz: 



D 








yu u uy 

■■'^2 C 











U^5i + s^52 cs{5i-52) 0^ 



cs{5i - ^2) s^5i + 0^82 
53 



(CI) 



(C2) 



(C3) 



where s (resp., c) denotes the sine (resp., the cosine) of angle 6 between axes x and X (see 
Fig.[I7D. 

Inserting the above equations into the evolution equation flTTl) . we obtain: 



c'/9i + s'(32 = 27cs(/3i - /^a) - + 3^2) 

^I(A^^ = 7(s'/5i + c^(32) - 2cs{5i - 82) 



(C4) 
(C5) 
(C6) 



Summing Eqs. flC4p and flC6P yields the evolution of Pi + ^2- The difference be- 
tween Eqs. (El) and ([06]), multiplied by (A - /32) cos(26'), plus Eq. ([CS]) multiplied by 
2(/3i — P2) sin(26'), yields the evolution of Pi — ^2- From there, the evolution of Pi and P2 
is readily obtained, see Eqs. (l66l) and ( !67j) . Using these equations, the difference between 
Eqs. (JCH) and then yields the evolution of cos (2^), see Eq. ([68]) 



2. Elastic shear flow 

Let the initial configuration be a material at rest, with Pi = P2 = ^- The system of 
Eqs. (l66l) . (l67j) and (168!) has a singularity at t = since Eq. (168!) contains a factor . In 
fact, the system of Eqs. flC4p . flC5l) and flC6P can be solved explicitely in the domain where 
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the material is purely elastic, with 6i = 82 = 0. Using the same notation as in Appendix I C It 

c2/5i + 5^/52 = I + 7V (C7) 
cs(A-/52) = 7t (C8) 
s^(3i + 0^(32 = 1 (C9) 



Whence: 



7(() = 7 * (CIO) 



ft(*) 



3. Plasticity onset 

When in the elastic regime, the system is described by Eqs.(...) above. Once the threshold 
is reached, given by Eq. (l82l) . the plastic term comes into play. 

Eq. (152]) by itself describes the limit of the elastic regime in terms of tensor B for any 
choice of the system history in the elasic regime (not just the simple shear implemented 
here) . It is represented on Fig. [191 



4. Stationary shear flow 

When the flow is stationary, the system of Eqs fl66H68|) can be simplified. In particular, 
Eq. (!68|) then implies: 

cos(2^) = u = (C14) 

Pi + P2 

=i„(29) = ^13^=^ (C15) 

Pi + P2 

and finally the tilt 6 of basis XY relative to basis xy can be expressed as: 

e = arctanj|^ (C16) 
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From Eqs (166|) and (167|) . still in the stationary regime, we obtain: 





(3i P2 



(C17) 



Combining this with Eq. (140|) which expresses the material incompressibility, we obtain 
53 = 0. 

Again from Eqs (!66ll and (!67l) in the stationary regime, we derive the shear rate: 



7 



7 



25i 
2(5i 



2^2 

>2 



2^2 



Using Eq. fIClSp . we then obtain: 



7 



(^1 - ^2 



(C18) 



(C19) 



Thus, the stationary branch of the flow plots are obtained by choosing pairs of values for 
/3i and (32 that satisfy Eq. (IC17p . where 5i and 82 are functions of (32 and (3^, = {^1^2)'^- 
To obtain the relevant values for /5i and (32, we first find the pair that verifies simulataneously 
the threshold condition given by Eq. (!82|) and the stationary condition given by Eq. (10170 . 
We then use a differential equation derived from Eq. flC17l) to follow the corresponding curve 
in the P1-P2 plane. 

Once the values for (3i and (32 are known, the shear rate is derived via Eq. ( lC19p . 



We can now slightly simplify Eqs. (184|) - (!9T|) in the stationary regime. 



a 



/3ia-l+/320'2 V/^ife (0-1-0-2) n 
/3i+/32 /3i+/32 
■yZ/^ife (0-1—02) /310-2+/320-1 g 



/3i+/32 









0"3 J 



(C20) 
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The first and second normal stress difference can be expressed as: 



Ni = ail - o"22 



No 



a ,n [«1 + «2(/3l + ^2)] 
Pi + P2 

^■22 - 0-3 

P1+P2 

1 



(C21) 



/9i + /32 /9i/52, 
1 



+^2 /3l/52 



(/5l/92) 



In the case of Mooney-Rivlin, coefficients ai and a2 are given by Eq. 
stress differences can be expressed as: 



N2 = ki 

+ /C2 



iPl - ^2? 
P1 + P2 

2/5i/?2 



ki + 



ko 



P1P2 



(C22) 

H) and the normal 
(C23) 



P1+P2 P1P2 



P1P2 P1+P2 



(C24) 



[1] L. D. Landau and E. M. Lifshitz. Theory of Elasticity. Butterworth-Heinemann, London, 3rd 
edition, 1995. 

[2] J. Friedel. Dislocations. Addison- Weskey Publishing Co., Lie, Reading MA, USA, 1983. 
[3] C. Kittel. Physique de I'etat solide. Dunod, Paris, 5th edition, 1983. 

[4] D. Weaire and S. Hutzler. The Physics of Foams. Oxford University Press, New York, 1999. 
[5] P. Coussot and C. Ancey. Rheophysique des pates et des suspensions. EDP Sciences, Les Ulis, 
Prance, 1999. 

[6] P. Coussot. Rheometry of pastes, suspensions and granular materials. Wiley, New York, USA, 
2005. 

[7] L. Cipeletti and L. Ramos. Slow dynamics in glasses, gels and foams. Curr. Opin. Colloid 
Interface Sci., 7:228-234, 2002. 



65 



[8] Michel Cloitre, Regis Borrega, and Ludwik Leibler. Rheological aging and rejuvenation in 

microgel pastes. Phys. Rev. Lett, 85(22) :4819-4822, 2000. 
[9] S. Cohen- Addad, R. Hohler, and Y. Khidas. Origin of the slow linear viscoelastic response of 

aqueous foams. Phys. Rev. Lett, 93:028302, 2004. 
[10] E. Eiser, F. Molino, G. Porte, and X. Pithon. Flow in micellar cubic crystals. Rheologica 

Acta, 39:201-208, 2000. 

[11] E. Eiser, F. Molino, G. Porte, and O. Diat. Nonhomogeneous textures and banded flow in a 

soft cubic phase under shear. Phys. Rev. E, 61:6759-6764, 2000. 
[12] J.-F. Berret, D. C. Roux, and G. Porte. Isotropic-to-nematic transition in wormlike micelles 

under shear. Europ. Phys. Journal E, 4:1261-1279, 1994. 
[13] J. P. Decruppe, S. Lerouge, and J.-F. Berret. Insight in shear banding under transient flow. 

Phys. Rev. E, 63:022501, 1999. 
[14] G. Porte, J.-F. Berret, and J. Harden. Inhomogeneous flows of complex fluids: Mechanical 

instability versus non-equilibrium phase transition. Europ. Phys. Journal E, 7:459-472, 1997. 
[15] C.-Y. David Lu, P. D. Olmsted, and R. C. Ball. Effects of nonlocal stress on the determination 

of shear banding flow. Phys. Rev. E, 84:642-645, 2000. 
[16] S. Lerouge, M. Argentina, and J. P. Decruppe. Interface instability in shear-banding flow. 

Phys. Rev. Lett, 96:088301, 2006. 
[17] L. Becu, S. Manneville, and A. Colin. Spatiotemporal dynamics of wormlike micelles under 

shear. Phys. Rev. Lett, 93:018301, 2004. 
[18] J.-F. Berret and Y. Serero. Evidence of shear-induced fluid fracture in telechelic polymer 

networks. Phys. Rev. Lett, 87:048303, 2001. 
[19] F. Molino, J. Appeh, M. Filah, E. Michel, G. Porte, S. Mora, and E. Sunyer. A transient 

network of telechelic polymers and microspheres: structure and rheology. J. Phys.: Condens. 

Matter, 12:A491-A498, 2000. 
[20] G. Debregeas, H. Tabuteau, and J.-M. di Meglio. Deformation and flow of a two-dimensional 

foam under continuous shear. Phys. Rev. Lett., 87:178305, 2001. 
[21] A. Kabla and G. Debregeas. Local stress relaxation and shear banding in a dry foam under 

shear. Phys. Rev. Lett, 90:258303, 2003. 
[22] M. L. Falk and J. S. Langer. Dynamics of viscoplastic deformation in amorphous solids. Phys. 

Rev. E, 57:7192-7205, 1998. 



[23] G. Picard, A. Ajdari, F. Lequeux, and L. Bocquet. Elastic consequences of a single plastic 
event : a step towards the microscopic modelling of the flow of yield stress fluids. Eur. Phys. 
J. E, 15:371, 2004. 

[24] G. Picard, A. Ajdari, F. Lequeux, and L. Bocquet. Slow flows of yield stress fluids: complex 

spatiotemporal behavior within a simple elastoplastic model. Phys. Rev. E, 71:010501, 2005. 
[25] M. Aubouy, Y. Jiang, J. A. Glazier, and F. Graner. A texture tensor to quantify deformations. 

Granular Matter, 5:67-70, 2003. 
[26] M. Asipauskas, M. Aubouy, J. A. Glazier, F. Graner, and Y. Jiang. A texture tensor to 

quantify deformations: the example of two-dimensional flowing foams. Granular Matter, 

5:71-74, 2003. 

[27] P. Marmottant and F. Graner. An elastic, plastic, viscous model for slow shear of a liquid 

foam. Eur. Phys. J. E, 23:337-347, 2007. 
[28] P. Saramito. A new constitutive equation for elastoviscoplastic fluid flows. J. Non-Newtonian 

Fluid Mechanics, 145:1-14, 2007. 
[29] O. Takeshi and K. Sekimoto. Internal stress in a model elastoplastic fluid. Phys. Rev. Lett, 

95:108301, 2005. 

[30] A. M. Kraynik and D. A. Reinelt. Microrheology of random polydisperse foam. In Proc. XIV^^ 

Int. Congr. on Rheology. The Korean Society of Rheology, 2004. 
[31] V. Mora. Etude de I'integration temporelle du tenseur taux de deformation. Application a 

la modelisation de I'elastoplasticite en grandes transformations. PhD thesis, Universite de 

Bretagne Sud, Lorient, FRANCE, 2004. 
[32] S. Benito et al. . in preparation. 

[33] R.J. Gordon and W.R. Schowalter. Trans. Soc. RheoL, 16:79-97, 1972. 

[34] A. Bertram. Elasticity and Plasticity of large deformations. Springer, Berlin, Heidelberg, New 
York, 2005. 

[35] R. Hohler and S. Cohen- Addad. Rheology of liquid foam. J. Phys.: Gondens. Matter, 
17:R1041-R1069, 2005. 

[36] R. Hohler, S. Cohen-Addad, and V. Labiausse. Constitutive equation to describe the nonlinear 
elastic response of aqueous foams and concentrated emulsions. J. RheoL, 48:679-690, 2004. 

[37] S. A. Khan, C. A. Schnepper, and R. C. Armstrong. Foam rheology: III. measurement of 
shear flow properties. J. RheoL, 32:69-92, 1988. 

67 



[38] F. Graner, B. Dollet, C. Raufaste, and P. Marmottant. Statistical tools to characterize discrete 
rearranging patterns, in 2 or 3 dimensions: cellular materials, assemblies of particles, preprint, 
2007. 

[39] J. Goyon, A. Colin ans G. Ovarlez, A. Ajdari, and L. Bocquet. Microfluidic velocimetry reveals 

spatial cooperativity in the flow of soft glassy materials, in preparation, 2007. 
[40] S. Marze, A. Saint-Jalmes, and D. Langevin. Protein and surfactant foams: linear rheology 

and dilatancy effects. Colloids and Surface A, 263:121-128, 2005. 
[41] J. Emile, E. Hardy, A. Saint-Jalmes, E. Terriac, and R. Delannay. Swelling of a single foam 

film under slipping. Colloids and Surface A: Physicochem. Eng. Aspects, 304:72-76, 2007. 
[42] E. Terriac, J. Etrillard, and I. Cantat. Viscous force exerted on a foam at a solid boundary: 

influence of the liquid fraction and of the bubble size. Europhys. Lett., 74:909-915, 2006. 
[43] C. R. Myers, B. E. Shaw, and J. S. Langer. Slip complexity in a crustal-plane model of an 

earthquake fault. Phys. Rev. Lett, 77:972-975, 1996. 
[44] Such interpretations were initially inspired by reflections on the dynamics of earthquakes 

(individual, local tectonic events relieve part of the stress locally but also report part of it to 



other places along the fault 



[45] In the presence of aging, randomly oriented Tl processes may sometimes locally enhance the 

stored deformation. On average, however, they are biaised by the ambiant stored deformation 

and thus tend to lower it ^. 
[46] As the additional stress that corresponds to the optional viscous element is transmitted to 

neighbouring regions of the material, it may induce additional stored deformation there, and 

hence impact the internal dynamics. 
[47] For real (positive or negative) values of the stress, Eq. ([T|) must contain an additional term: 

G 

& = G D - — {{a - ay) 6[a - ay] + {a + ay) 9[-a - ay]) 

[48] The deformation e defined by Eq. (jl5|) happens to coincide, in the present case where F is 
symmetric (see Fig. [3]), with the usual (Lagrangian) deformation (as defined for instance by 
Landau [l|), relative to the relaxed state: 

1 



^Landau 



-{F^F-l) 



[49] This decomposition results from the Theorem of representation of isotropic functions. 
[50] up to the effect of disorder, see paragraph IVI C[ 

68 



[51] Note that this elastic compressibility can always be defined on short time scales, even if the 

threshold stress for plasticity is arbitrarily small. 
[52] We assume that all non-linearities in D are included in the dependence on the deformation 

B, which itself results from the applied deformation rate D. 
[53] See Appendix IC 21 for the apparent singularity at t = when Pi = /?2 and for the state of the 

system in the elastic regime. 
[54] a formulation specific to 2D materials. 

[55] reduction of the stored deformation while the total deformation rate is zero. 
[56] To derive Eq. (jl7p from Eq. ()A2p . we use the identity: 

= d[B- B-^]/dt = B • d[B~^]/dt + [dB/dt] • B-^ 



69 



